function [Ss, rs] = find_spillovers_dyn(Ss, lambdas, rs, Prob, tol,max_It)
    global bet0 bet1 tau a b w e pi_theta N_a N_theta N_w N_e w_max w_min alpha myslope...
        xi xi1 ub lb pi_e H w theta Int_pi_w sigmae elasts Pr log_e_mu log_e_sigma phi sigma gamma rts r_bar fp r_base w_mult w_bar...
        
   
    % Fixed Point Iteration

    pi_w = sum(Prob, 2);
    w_bar = dot(w, pi_w ); 

        
    for it = 1:max_It
        % Spectral Algorithm to find market clearing prices
        x_0 = rs(1:end);
        x_1 = x_0 + (size(rs,1):-1:1)'/100; 

        f_0 = find_rents_dyn(x_0, Ss,  lambdas, Prob);
        f_1 = find_rents_dyn(x_1, Ss,  lambdas, Prob);

        max_iter = 100;
        for i = 1:max_iter 
            s = x_1(1:end) - x_0(1:end);
            y = f_1(1:end) - f_0(1:end);

            x_1 = x_0;
            f_1 = f_0;

            lambda = dot(s,y)/dot(y,y);
            
            step_size = 1;
            for iter_sub = 1:20
                
                x_0_t = x_1 - lambda * f_1 * step_size;
          
                if sum(x_0_t(1:end-1) > 0 ) == (length(x_0_t)-1)
                    break
                else
                    step_size = step_size * 0.9;
                end
            end
            x_0 = x_0_t;
            f_0 = find_rents_dyn(x_0, Ss,  lambdas, Prob);

           
            diff = sqrt(dot(x_0 - x_1, x_0 - x_1));
            diff_f = sum(abs(f_0));
            
            if (diff < 1e-6) & (diff_f < 1e-3)
                rs = x_0;
                break
            elseif i == max_iter
       
                rs = [w_bar*0.7, w_bar*0.5, w_bar*0.2]';
                x_0 = rs(1:end);
                minf = @(x) sum((find_rents_dyn(x, Ss,  lambdas, Prob)).^2);
                options = optimset('TolFun',1e-6,'TolX',1e-3);
                [minx, fval, exitflag] = fminsearch(minf, x_0(1:end), options );
                
                rs = minx; 
                if (fval > 1e-3 | isnan(fval) | (exitflag ~= 1)) & (it > 5 )
                    fval, minx
                    error("Did not find root")
                end
            end                
        end
      
        % Calculate Choice Probabilities
        [P_S, ~, wprimeS, ~] = choiceProb(Ss, rs);
        P_S = bsxfun(@times, P_S, reshape(pi_theta, [1,1,1,2]));
        popMass = bsxfun(@times, P_S, Prob);
        
        wMass = squeeze(sum(sum(bsxfun(@times, wprimeS, popMass), 1:2), 4));
        popMass = squeeze(sum(sum(popMass, 1:2),4));
        
        Ssi = wMass ./ popMass;       
        Ssi = Ssi.^xi1;
        
        d_S = abs(Ssi - Ss);

        % Test for convergence
        if sum(d_S) < 2 * tol
            break
        elseif it == max_It
            display(sum(d_S(:)))
            error("Did not converge in Spillover") 
        end
        Ss = Ssi;
     
    end
    
end